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Wolbachia causes the feminization of chromosomally male embryos in several species 
of crustaceans and insects, including the leafhopper Zyginidia pullula. In contrast 
to the relatively well-established ecological aspects of male feminization (e.g., sex 
ratio distortion and its consequences), the underlying molecular mechanisms remain 
understudied and unclear. We embarked on an exploratory study to investigate the 
extent and nature of Wolbachia's effect on gene expression pattern in Z pullula. 
We sequenced whole transcriptomes from Wolbachia-Mected and uninfected adults. 
18147 loci were assembled de novo, including homologs of several Drosophila sex 
determination genes. A number of transcripts were flagged as candidate Wolbachia 
sequences. Despite the resemblance of Wolbach /a-infected chromosomal males to 
uninfected and infected chromosomal females in terms of sexual morphology and 
behavior, principal component analysis revealed that gene expression patterns did not 
follow these sexual phenotype categories. The principal components generated by 
differentially expressed genes specified a strong sex-independent Wolbachia effect, 
followed by a weaker Wolbach ia-sexual karyotype interaction effect. Approaches to further 
examine the molecular mechanism of l/l/o/tec/7/a-host interactions have been suggested 
based on the presented findings. 



Keywords: Wolbachia infection, male feminization, principal component analysis (PCA), Zyginidia pullula 
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INTRODUCTION 

Wolbachia is an intracellular symbiont alpha-proteobacterium 
that infects a wide range of arthropods and nematodes 
(Schulenburg et al, 2000; Werren et al., 2008). It is often trans- 
mitted vertically from females through the eggs to their future 
progeny; although, horizontal transfer between hosts has also 
been documented (Werren et al, 1995; Cordaux et al., 2001). 
Studying the mechanism of Wolbachia-host interactions is fasci- 
nating for many reasons. Wolbachia is capable of inducing several 
intriguing sex-related phenotypes in its hosts, including male 
killing (MK), in which infected males die during embryonic or 
larval stages; male feminization (MF), that is the development 
of genetic males into females; thelytokous parthenogenesis (TP) 
in which infected virgin females produce daughters. All of these 
phenotypes distort the progeny sex ratio in favor of females thus 
ensuring higher transmission rate of Wolbachia to the next gen- 
eration of hosts (Werren et al., 2008; White et al, 2013). Another 
fascinating effect of the infection is cytoplasmic incompatibility 
between gametes (CI), which results in aberrant or consider- 
ably reduced offspring production, if uninfected females mate 
with infected males, or if the parents are infected with differ- 
ent Wolbachia strains (Werren et al, 2008; White et al, 2013). In 
this case, infected females possess a reproductive advantage com- 
pared to uninfected ones, and this again ensures the spreading of 



Wolbachia into the host population. Fast transition between the 
four phenotypes in the course of the coevolution of Wolbachia 
and its hosts hints that similar molecular mechanisms might 
underlie the apparently different effects (Ma et al., 2014). Due to 
its enormous host range, Wolbachia may have played a crucial role 
in the evolution of sex determination system and reproductive 
strategies in arthropods (Cordaux et al., 2011; Awrahman et al., 
2014; Ma etal, 2014). 

Various approaches have been employed to investigate the 
Wolbachia-host interactions in naturally infected and uninfected 
strains (Hoffmann et al, 1990; Negri et al, 2006; Riparbelli et al., 
2012), experimentally inoculated cell lines (Noda et al., 2002; Xi 
et al., 2008), and antibiotic treated specimens (Hoffmann et al., 
1990; Casiraghi et al., 2002). Although Wolbachia is an obligate 
intracellular symbiont natuarally, protocols have been developed 
to keep it viable in cell-free media for days; however, no repli- 
cation occurs in the extracellular phase (Rasgon et al., 2006; 
Gamston and Rasgon, 2007). The experimental/analytical tech- 
niques comprised a wide range including classical crossing and 
fecundity measurements (e.g., Hoffmann et al, 1990; Dunn et al., 
2006), microscopic approaches (in situ hybridizations, electron 
microscope and immunohistochemical techniques for bacterium 
detection inside hosts and cells, tissues, etc.) (e.g., Negri et al., 
2008; Fischer et al., 201 1), gene expression analysis (e.g., Xi et al., 
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2008; Kremer et al., 2009, 2012; Hughes et al, 2011; Chevalier 
et al., 2012; Darby et al, 2012; Liu et al, 2014), bioinformatic 
genome sequence annotation and functional prediction (e.g., Wu 
et al., 2004; Foster et al, 2005; Klasson et al., 2008), and math- 
ematical modeling of the ecological consequences of CI or sex 
ratio distortion (e.g., Taylor, 1990; Turelli, 1994). Despite all these 
efforts, a coherent mechanistic story of Wolbachia's effect is still 
lacking. The picture is incomplete even for CI which occurs 
in Drosophila and is the most extensively studied Wolbachia- 
induced phenotype; although, cytoskeleton reorganization and 
asynchrony in nuclear envelope break down and chromosomal 
condensation of male and female pronuclei after fertilization 
have been implicated in the process (Serbus et al., 2008; Werren 
et al., 2008). The other three phenomena are less well under- 
stood. TP seems to result from induction of diploidy in species 
with a haplodiploid sex determination system by production and 
development of diploid eggs; that is achieved by altering meio- 
sis to produce diploid gametes (Weeks and Breeuwer, 2001), the 
abortion of the first mitotic division after chromosomal dupli- 
cation (Pannebakker et al., 2004), or the fusion of the two 
haploid nuclei after first mitosis of induced eggs (Gottlieb et al., 
2002). The molecular bases of MK and MF are least understood 
but they are suspected to share certain components as MK is 
often the result of a lethal and incomplete attempt at feminiza- 
tion of genetic male embryos (Werren et al, 2008). The most 
direct mechanistic evidence comes from the study of male killing 
Wolbachia in the moth Ostrinia scapulalis showing that it over- 
rides the karyotypic signal in genetic males to produce the female 
dsx isoform (Sugimoto and Ishikawa, 2012). This suggests that 
Wolbachia impacts the sex determination pathway at or above 
dsx. Apart from this direct effect on the pivotal sex determin- 
ing gene dsx, MK or MF Wolbachia infection is reported to be 
accompanied with defective chromatin remodeling (Riparbelli 
et al., 2012), induction of host immune response (Chevalier et al., 
2012), and epigenetic reprogramming of the host (Negri et al., 
2009a). 

Zyginidia pullula is a leafhopper with XX/XO male heteroga- 
metic sex determination system in which Wolbachia causes fem- 
inization of chromosomal males (Negri et al, 2006). Infected 
female leafhoppers are morphologically indistinguishable from 
uninfected females; but feminized chromosomal males have an 
intersex phenotype i.e., they have the upper pygofer appendages, 
a typical male secondary sexual feature. These appendages show 
varying degrees of development, from being fully developed in 
some specimens to being a barely recognizable stump in oth- 
ers (Negri et al, 2006). Feminized males with upper pygofer 
appendages reduced to a stump have ovaries morphologically 
similar to uninfected females, whereas those with prominent 
appendages possess malformed and probably less functional 
ovaries (Negri et al., 2008). The "degree of feminization" has 
been shown to be correlated with Wolbachia density in the host 
tissues in several systems (Jaenike, 2009). We have previously 
reported that Wolbachia instigates epigenetic reprogramming of 
Z. pullula (Negri et al, 2009a,b) and probably interacts with the 
insect hormone biosynthesis pathway to stimulate the produc- 
tion of feminizing hormones (Negri et al., 2010; Negri, 2012). 
In this study, whole transcriptomes of male and female Zygindia 



samples ( Wolbachia-infected and uninfected) were analyzed with 
Illumina deep sequencing technique, in order to understand 
the scope and nature of the Wolbachia-induced change in the 
host gene expression profile. Our initial idea was that if male 
feminization is the main consequence of Wolbachia infection, 
transcriptomes from the three female types (uninfected females, 
infected females and feminized males) should resemble each 
other and be different from the only phenotypically male group 
(uninfected males). In fact, we decided to test the hypothesis 
that sex reversal is Wolbachias main effect at the transcriptome 
level. Were this confirmed, we would proceed to identify dif- 
ferentially expressed genes between the two sexual phenotype 
groups. 

METHODS 
ZYGINIDIA SPECIMENS 

34 overwintering females of Z. pullula were collected in the same 
grass field in north Italy; and were reared individually in the lab- 
oratory as described in Negri et al. (2006). Overwintering females 
have often mated with several males (rarely with only one). 
By carefully examining the progeny, Wolbachia-infected (i.e., all 
female brood) and uninfected (i.e., male and female brood) lines 
were identified. Wolbachia infection was then confirmed by PCR 
on the mothers and randomly chosen samples from the brood 
as described in Negri et al. (2006). Morphological investigation 
as to the presence or absence of upper pygofer appendages lead 
us to separate feminized males from genetic females in the all- 
female (i.e., Wolbachia-infected) lines, and males and females in 
the uninfected lines. Males from uninfected lines were mated to 
the physiologically female progeny of the infected lines (consist- 
ing of genetic females and males) at each generation to produce 
the next generation of infected females (and feminized males). 
This backcrossing to uninfected males was done for at least three 
generations in the lab. Fifty adults from each of the four differ- 
ent categories of uninfected females (F), uninfected males (M), 
infected females (FW) and feminized (infected) males (MW) were 
pooled together for RNA sequencing. 

cDNA LIBRARY PREPARATION AND SHORT-READ SEQUENCING 

cDNA libraries were made from male and female specimens 
of infected and uninfected leafhopper lines. Infected males are 
phenotypically intersex and exhibit different degrees of feminiza- 
tion depending on the concentration of Wolbachia, ranging from 
individuals with functional ovaries to individuals with female 
secondary sexual characters, but possessing testes. We used thor- 
oughly feminized infected males for RNA extraction. RNA purifi- 
cation, cDNA synthesis and Illumina library construction were 
performed using the protocols of Mortazavi et al. (2008), with the 
following modifications: total RNA, mRNA and DNA were quan- 
tified using a Qubit fiuorometer (Invitrogen); mRNA fragmenta- 
tion was performed using Fragmentation Reagent (Ambion) for 
a 3 min and 50 s incubation at 70° C and subsequently cleaned 
through an RNA cleanup kit (Zymo Research); additional DNA 
and gel purification steps were conducted using Clean and 
Concentrator kits (Zymo Research). Each sample library was 
sequenced as pair-ended 76-base reads on an Illumina Genome 
Analyzer II. 
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DE NOVO TRANSCRIPTOME ASSEMBLY AND EXPRESSION LEVEL 
CALCULATION 

Due to the sensitive nature of de novo assembly, it is critical that 
the reads used to generate contigs have the highest sequencing 
quality. Reads were removed from consideration in the de novo 
assembly if they had a terminal phred (Ewing and Green, 1998) 
quality value less than 15, or contained more than 2 unknown 
nucleotides (i.e., N). Reads were also filtered due to similarities 
to known PCR primer and Illumina Adapter sequences. Using 
the reads pooled from all of the four samples that were not fil- 
tered out, the de novo assembly program Velvet (version 1.0.15) 
(Zerbino and Birney, 2008) was used in conjunction with a cus- 
tom post-processing algorithm capable of retaining information 
from alternative splices (Sze et al., 2012) to assemble short reads 
into contigs, using sequence overlap information until the con- 
tigs could no longer be extended. Velvet was run under the 
following settings with a kmer length of 35: -cov_cutoff auto 
-max_branch_length o -max_divergence 0 -max_gap_count 0 - 
read_trkg yes. Sequenced reads that were kept as pairs and not 
filtered out together or separately were treated as "-shortPaired" 
with insert length of 175 bases and standard deviation of 75 bases. 
Single end reads that were not filtered out were treated as "-short." 

With the set of de novo assembled sequences serving as 
a reference, reads from each of the individual samples were 
mapped using the Burrows-Wheeler Aligner (BWA) (Li and 
Durbin, 2009). The number of reads that mapped to the con- 
tigs of each gene was tabulated and normalized to calculate 
FPKM (Fragments Per Kilobase Of Exon Per Million Fragments 
Mapped). Additional normalization among all samples was per- 
formed using the TMM protocol (Trimmed Mean of M-values) 
outlined in Robinson and Oshlack (2010), which takes into 
account differences in overall RNA populations across samples 
and is one of several methods used to evaluate RNA sequencing 
data. Normalization was implemented using the edgeR package in 
R (Robinson et al, 2010). All statistical analyses and graphs evalu- 
ating consistency between samples were produced using R v2. 1 3.0 
(R Development Core Team, 201 1). 

GENE FUNCTIONAL ANNOTATION AND CLASSIFICATION 

Blast2GO v.2 (Gotz et al., 2008) and WEGO (Ye et al, 2006) 
were used to obtain Gene Ontology (GO) annotations. Genes 
were also annotated using a BLASTX search (Altschul et al, 1990) 
(Expected value <1.00e-05) to the nr protein database available 
from GenBank as well as to the set of protein sequences avail- 
able from the Drosophila melanogaster 5.34 and the pea aphid 
Acyrthosiphon pisum 2.1 releases. We chose the annotation with 
the highest BLAST score as long as the span of the alignment was 
greater than 80% of the length of the contig under query. For 
genes that did not report any hits, we lowered the minimum span 
to 40% of the length, choosing the annotation with the highest 
BLAST score having Expected value < 1.00e-05. 

PRINCIPAL COMPONENT ANALYSIS OF GENE EXPRESSION VALUES 

Expression values were cleaned of extreme outliers, quartile- 
normalized and log-transformed before they were used for PCA. 
To make sure the result were not artifacts of the data prepara- 
tion method, PCA was repeated on the raw (not normalized, not 



log-transformed) expression values as well as after several differ- 
ent outlier-filtering and normalization strategies. These statistical 
procedures were done in SAS 9.3. 

RESULTS 

SHORT-READ SEQUENCING AND DE NOVO ASSEMBLY 

The mRNA population was analyzed with Illumina deep sequenc- 
ing of male and female Zygindia samples with and without 
Wolbachia infection. The pooled data from all samples had a total 
of 50 M pair-ended reads that were 76 bases long. All Illumina 
sequences are available for download at the NCBI Short Read 
Archive under the BioProject PRJNA171390. After sequences were 
filtered based on quality and matches to adapter and primer 
sequences, the 38 M reads from all four samples were pooled 
together and run through Velvet and the post-processing algo- 
rithm. Eventually, 18,147 loci and a total of 27,236 transcripts 
were assembled; multiple transcripts of a locus pertained often 
to different splicing isoforms and occasionally to largely differ- 
entiated alleles. The transcripts ranged in lengths from 291 bp 
to 15,389 bp, with mean and median lengths of 1006 bp and 
702 bp, respectively. This assembly included a fairly large num- 
ber of long transcripts: 25% were longer than 1250 bp and 10% 
were longer than 2000 bp. Of the 18,147 loci, 14,068 (77.5%) had 
a single isoform and the remaining 22.5% had multiple ones. 
Transcripts within a locus were subsequently collapsed into a sin- 
gle "representative locus sequence" by using ClustalW to run a 
multiple sequence alignment and identifying the locus consen- 
sus sequence. Mean and median lengths of consensus sequences 
were 900 bp and 618 bp, respectively. The total length of all loci 
consensus sequences was 16.3 Mb. 

GENE FUNCTIONAL ANNOTATION AND CLASSIFICATION 

6946 loci, corresponding to 38% of the entire dataset, were Gene 
Ontology annotated with Blast2GO. The consensus sequences 
were also aligned using a BLASTX search to the nr protein 
database available from GenBank as well as to the set of protein 
sequences available from Flybase and the aphid genome. Table 1 
shows the proportion of cases that resulted in a hit where the 
length of the alignment was greater than 80% or 40% of the length 
of the query (leafhopper sequence). One might very crudely 
attribute the 80% alignment span hits to true genie homology and 
the 40% alignment span hits to conserved domains. 



Table 1 | Summary statistics of Zyginidia transcripts homology 


search. 








40% homology length 


80% homology length 


Genbank 


60% 


39% 


Flybase 


48% 


26% 


Pea aphid 


81% 


32% 



40% homology length indicates that the length of the homologous segment cov- 
ers at least 40% of the query (leafhopper) sequence. A correspondingly similar 
definition applies to the 80% homology length category. The percent values in 
the table cells show what percent out of all loci (18147) fit each criterion when 
Blasted against the designated dataset. 
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A number of genes potentially involved in the leafhopper sex 
determination were identified through homology search with 
the Drosophila sex determination genes. Although pea aphid is 
Zyginidia s closest relative with a reference genome sequence (The 
International Aphid Genomics Consortium, 2010), the func- 
tional annotation for this genome is not as complete as that of 
Drosophila. Sex determining genes of pea aphid have been found 
based on homology with Drosophila sequences and lack direct 
experimental verification (The International Aphid Genomics 
Consortium, 2010). Therefore, we decided to use Drosophila 
sequences as the reference set. Figure 1 depicts the canonical 
sex determination pathway in Drosophila. Homologs of several 
Drosophila sex determination genes were identified among the 
transcripts including dsx (doublesex), tra-2 (transformer- 2), vir 
(virilizer), fl(2)d (female lethal d), snf (sans fille) and ix (inter- 
sex). No leafhopper homologs could be identified for tra (trans- 
former), sxl (sex lethal), fru (fruitless) or her (hermaphrodite). 
Table 2 shows the expression levels for the identified leafhopper 
sex determination genes. 

Seventeen genes in our dataset were flagged as likely Wolbachia 
sequences according to the Blast results against the NCBI dataset. 
Bacterial origin seems very probable for a number of these 
transcripts based on the expression levels in infected and unin- 
fected lines, plus high similarity to known Wolbachia sequences 
(Table 3). These sequences were Blasted against the aphid genome 



Femaleness 

XX; AA 
Vir, FL(2)D, Snf j 

I 




Female CNS Female soma 

X / 

Female sexual 
behaviour 

FIGURE 1 | Sex determination pathway in Drosophila, modified from 

Sanchez (2008). Sxl F and Sxl M refer to functional female and nonfunctiona 
male isoforms of the Sxl protein, respectively. Tra F is the functional female 
form of the Transformer protein, which in conjunction with the constitutive 
gene product Tra-2 controls female-specific splicing of dsx and fru. snf, vir 
and FL(2)D are required for late female-specific splicing of Sxl but play no part 



to check if there was an indication of horizontal transfer; they 
were also Blasted against the Drosophila genome as a distant 
outgroup (Table 3). 

PRINCIPAL COMPONENT ANALYSIS OF TRANSCRIPTOMES 

Principal component analysis on transcriptomes of the four 
leafhopper samples surprisingly revealed that Wolbachia infec- 
tion changes the host transcriptome extensively and the effect 
is by no means limited to sex-reversal. As evident in Figure 2, 



Table 2 | Homologs of Drosophila sex determination genes in the 
Zyginidia transcriptome and their normalized expression levels. 



Locus 


Fly homolog 


Multiple 
isoforms? 


FW 


MW 


F 


M 


5652 


dsx 


Y 


11.6 


49.4 


11.7 


19.7 


8743 


tra-2 


Y 


5.9 


13.2 


8.8 


26.1 


5015 


vir 


N 


20.2 


9.9 


16.1 


23.2 


10229 


fl(2)d 


N 


17.4 


32.9 


40.9 


63.3 


21060 


snf 


Y 


0 


3.3 


20.5 


10.9 


18743 


ix 


N 


5.9 


11.6 


81.8 


43.9 



FW, infected female; MW, infected feminized male (intersex female); F, unin- 
fected female; M, uninfected male. 
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in determining early Sxl splicing pattern. The genes for which Z pullula 
homologs have been identified in this study, are boxed in gray. For more 
details on the regulation and function of these genes, refer to (Sanchez, 
2008; Gempe and Beye, 201 1 ) or other similar resources. Reproduced with 
permission from The International Journal of Developmental Biology (Int. J. 
Dev. Biol.) (2008) Vol:52, pp. 837-856. 
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the first PC (explaining 66.46% of variance) is highly corre- 
lated with all of the samples indicating that the expression of 
most genes is not significantly altered by Wolbachia and is sim- 
ilar across all samples. By the second PC (explaining 20.36% of 
variance), Wolbachia infected male and female samples cluster 
together and uninfected male and female cluster together. This PC 
is generated by genes whose expression is changed by Wolbachia 
consistently regardless of sexual karyotype or phenotype. The 
third PC (explaining 7.97% of variance) indicates an interaction 
term: F and M are similar and stand in the middle of the scale, 
with MW and FW occupying the opposite sides of them. This PC 
is generated by genes that are expressed similarly in uninfected 
males and females, and Wolbachia infection changes their expres- 
sion in opposite ways in chromosomal males and females. Overall, 
sex inversion does not seem to be the only or even the biggest 
effect of Wolbachia on gene expression patterns in Zyginidia, even 
if it is the most conspicuous phenotypic consequence; otherwise, 
we would expect the three phenotypically female groups (F, FW 
and MW) to cluster together and the only male group (M) to 
stand separate from them. None of the PCs show such a pattern. 
PCA was repeated on expression values without the initial outlier 
filtering, and applying several different normalization and trans- 
formation strategies; they all yielded the same picture as described 
above: the main effect was invariably the presence or absence of 
Wolbachia regardless of sex (details not shown). 

DISCUSSION 

We assembled the Z. pullula transcriptome de novo and pro- 
duced 18,147 loci and 27,236 transcripts with a total consensus 
sequence length of 16.3 Mb. These numbers were well within 
the expected range based on the aphid genome information. The 
aphid genome was reported to contain 11,089 highly supported 
RefSeq gene models with a total exonic length of 21.6 Mb; adding 
the gene models from six other gene prediction programs, a total 
of 34,604 non-redundant gene models with the total exonic length 
of 35.7 Mb were described (The International Aphid Genomics 
Consortium, 2010). The true number of genes is purportedly a 
number between those two estimates. Hence, our de novo assem- 
bly of the transcriptome seems to have captured a reasonable 
proportion of the expressed genes. 

The results of sequence homology search (Table 1 ) confirm the 
closer relatedness of Z. pullula to A. pisum (the aphid) than to 
Drosophila. A caveat to this analysis is the extensive set of dupli- 
cations in the aphid genome (The International Aphid Genomics 
Consortium, 2010). Without a leafhopper reference genome, we 
do not know if the same wave of duplications has affected Z. pul- 
lula or not; however, there was an indication in our data that it 
might have. By visual inspection of the sequences that were anno- 
tated as isoforms of a single locus computationally, we realized 
that some of them did not show signatures of known alterna- 
tive splicing patterns; but looked like highly differentiated alleles 
(details not shown). These may indeed be paralogous sequences 
in the process of divergence. Further investigation, including the 
sequencing of single individuals rather than pools of them, will be 
required to separate paralogy from allelic variation. 

A number of leafhopper sex determination genes were 
identified based on homology with fly sequences (Table 2). 
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FIGURE 2 | Principal component analysis of gene expression levels in 
four Zyginidia samples: infected females (FW), infected males or 
intersex females (MW), uninfected females (F) and uninfected males (M). 

The first three PCs are depicted here. PC1 is correlated with all sample labels 
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but PC2 separates the infected and uninfected samples conclusively. PC3 
reveals an interaction term between Wolbachia status and sexual karyotype. 
The "_ql" suffix after line names means that the expression values were 
quantile-normalized and log transformed. (A) PC1 vs. PC2; (B) PC1 vs. PC3. 



Insect sex determination machinery has evolved around the 
tmnsformer-doublesex axis (Sanchez, 2008); tra is the fast evolving 
component responsible for receiving the signal-sometimes 
through mediators- from the upstream sex determining factors 
(chromosomal constituent, incubation temperature, etc.), and 
dsx is the conserved switch that relays this signal down to the 
developmental processes (Sanchez, 2008; Verhulst et al., 2010). 
It is, therefore, not surprising that we found a homolog for dsx 
and not for tra in our dataset. The short length of the aligned 
segments prevented reliable assignment of male and female iso- 
forms; but these initial results can be used to design primers 
to extract the whole genes from the leafhopper genome. Future 
experiments can then follow the flow of the signal in the sex deter- 
mination pathway to identify where the cascade is diverted to 
female development in Wolbachia-infected genetic males. In the 
moth O. scapulalis, the impact point is somewhere above the level 
of dsx (Sugimoto and Ishikawa, 2012). Having the sequences of 
dsx male and female isoforms, one could check whether this is also 
true in leafhoppers. Unfortunately, the lack of replicates in our 
preliminary data makes it impossible to assess the significance of 
differential expression of genes across our four groups (FW, MW, 
F, and M). This is another task that remains to be done in future 
projects. In addition, development of X-linked sequence markers 
will enable early sexing of the embryos (based on the female XX / 
male XO karyotypes) through quantitative PCR; and facilitate the 
study of early developmental processes in infected and uninfected 
specimens. 

We found a number of Wolbachia-related transcripts in the 
sequenced cDNA libraries (Table 3). The loci expressed mainly in 
infected lines with great similarity to known Wolbachia sequences 
are likely to have Wolbachia origin (e.g., loci 1053, 1097, 1331, 
and 13961). Curiously, a couple of loci are expressed primarily 



in the uninfected lines (e.g., locus 22635). At this point, we 
do not have a hypothesis as to the reason behind this obser- 
vation. Repeating the experiments with replicates and higher 
sequencing depth would be the first step to confirm the repro- 
ducibility of these patterns. Our protocol of mRNA purification 
for creation of cDNA libraries involved a hybridization step with 
oligo-T ligands, which targets the eukaryotic mRNA poly-A tails; 
therefore, it will be necessary to employ a different purifica- 
tion strategy in order to capture most of the poly-A lacking 
bacterial mRNAs. Table 3 shows that several of the Wolbachia- 
related sequences code for Ankyrin-repeat proteins. Wolbachia 
genomes are well known for containing an extraordinarily high 
number of these genes (Wu et al., 2004; Iturbe-Ormaetxe et al., 
2005). Gene transfer between Wolbachia and mosquito hosts 
has been previously reported (Woolfit et al, 2009). PCR exper- 
iments and phylogenetic analyses have confirmed horizontal gene 
transfer from bacterial endosymbionts to the aphid genome 
(The International Aphid Genomics Consortium, 2010). Similar 
approaches will be required to confirm bacterial or insect ori- 
gin for the transcripts listed in Table 3. We tried to check 
for possible aphid lineage-specific horizontal transfers by ask- 
ing whether a likely Wolbachia transcript shows high sequence 
similarity to an aphid sequence, but not a fly sequence; none 
of the loci in Table 3 expressed such a pattern. One of the 
Wolbachia-related transcripts showed a degree of homology with 
the aphid vasa gene (locus 4382). Almost identical homologs of 
this sequence exist in the three published Wolbachia genomes 
(Blast results not shown); its homologs in fly, leafhopper and 
the published Wolbachia genomes are characterized or predicted 
ATP-dependent RNA helicases. vasa has been implicated in trans- 
mission of maternal effects and sex determination in clams 
(Milani et al, 201 1). It will be very interesting to check if products 
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of host-homologous genes are actually exported out by Wolbachia 
into the host cell. 

We used natural isolates of infected and uninfected leafhoppers 
for our comparisons with no antibiotic treatment. This relieved 
our comparisons from the confounding effects of antibiotic treat- 
ments on the host physiology. The rationale behind the tradi- 
tional use of antibiotics to cure the infected lines from Wolbachia 
is to obtain infected and uninfected lines with the same genetic 
background. However, antibiotics can change the host physiology 
substantially, and quite remarkably, their effect can perpetuate 
through several generations of unexposed progeny (Ballard and 
Melvin, 2007; Zeh et al., 2012; Fridmann-Sirkis et al, 2014). We 
avoided the use of antibiotics completely and achieved homoge- 
nous genetic backgrounds among samples by taking advantage of 
repeated backcrossing of infected females to uninfected males. We 
collected all of our founder specimens from the same leafhop- 
per population in a grass field. In the sampled population, the 
sex-ratio was only moderately female biased, with a moderate 
prevalence of the infection ( — 1:1.8 male:female, Wolbachia infec- 
tion rate ~30% of the collected females; Negri I., unpublished 
data). As uninfected males are the only physiological males in 
existence, all the "egg-laying females" (in the field and in the lab, 
including the females used in this study) always mate with (and 
only with) uninfected males. Thus, all of our infected and unin- 
fected lines come from the same genetic background. We carried 
out three further generations of backcrossing of infected females 
to uninfected males in the lab to effectively remove any residual 
genetic variation between the two groups. Details of rearing con- 
ditions are described in Negri et al. (2006). The natural pattern of 
sexual reproduction and the additional backcrossing done in the 
lab ensure the similarity of nuclear genetic backgrounds. We also 
tested mitochondrial gene sequences in Zyginidia samples from 
different Italian localities, both infected and uninfected, and they 
were all nearly identical (Negri I., unpublished data). 

Through principal component analysis, we have showed that 
Wolbachia-'mduced changes in the host transcriptome are mainly 
sex- independent, and cannot be explained only by the sex reversal 
of genetic males. Previous transcriptomic studies on Wolbachia 
have reported changes in the expression of genes unrelated to 
the reproductive phenotype. For instance, Wolbachia infection in 
Armadillidium vulgare triggered the overexpression of immune- 
related genes (Chevalier et al., 2012). In the parasitoid wasp 
Asobara tabida, endosymbiont infection or lack thereof was asso- 
ciated with changes in expression of genes related to female repro- 
ductive development, iron and oxidative stress regulation, and 
immune recognition (Kremer et al., 2009, 2012). Artificial infec- 
tion of Anopheles cell cultures by Wolbachia, surprisingly caused 
down-regulation of immune, stress response and detoxification 
genes (Hughes et al., 2011). Wolbachia-moculated Drosophila cell 
lines exhibited differential expression of several GO categories not 
directly related to reproduction, including antimicrobial humoral 
response, ion homeostasis, response to unfolded protein and 
response to chemical stimulus (Xi et al, 2008). In Aedes aegypti, 
Wolbachia was shown to manipulate the expression of a met- 
alloprotease gene through induction of a specific host miRNA 
(Hussain et al., 201 1). Apart from such direct evidence, the obser- 
vation of various forms of fitness cost in the feminized males, is 



consistent with the idea that sex reversal is not the sole effect of 
feminizing Wolbachia (Moreau et al., 2001; Rigaud and Moreau, 
2004). Nevertheless, our study is the first one to quantitatively 
demonstrate that infection itself has a larger effect than that of 
sex reversion, through PCA of all of the available gene expression 
levels. 

Lack of replicates meant that we could not quantitatively iden- 
tify differentially expressed genes between the lines because we 
could not calculate variances. Instead, we focused on the global 
patterns of gene expression by applying PCA to gene expression 
values. Thousands of loci (each acting as one observation point) 
were used to generate the PCs. Antibiotic treatment and differ- 
ent genetic backgrounds could have been two potential sources of 
systematic bias in this type of analysis; they could have generated 
similar clustering patterns and confounded the interpretation of 
results. However, through the single-population sampling and 
the repeated backcrossing scheme, we avoided both sources of 
confusion. 

Based on the PCA results, we encourage the use of biochemi- 
cal bottom-up approaches focusing on the whole Wolbachia effect 
rather than the specific sex inversion event. Wolbachia's effect 
is perceivably mediated by molecules secreted into the host cell 
or expressed on the outer membrane surface of the bacterium- 
containing vesicles. Wolbachia cannot be maintained in cell-free 
cultures indefinitely; but there are protocols to keep them alive in 
synthetic media for several hours (Rasgon et al., 2006; Gamston 
and Rasgon, 2007). In such a setting, the molecules released into 
the medium can be detected and purified using chromatographic 
and/or mass spectrometric approaches. Appropriate methods 
can be used, too, for isolation and characterization of surface 
molecules from the bacterium-containing vesicles. Pull-down 
experiments on the host proteins by these Wolbachia released or 
surface molecules might reveal the initial cellular targets of the 
endosymbiont-host interaction. 
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